Differential regulation of mRNAs and lncRNAs related to lipid metabolism in Duolang and Small Tail Han sheep

The function of long non-coding RNA (lncRNA) can be achieved through the regulation of target genes, and the deposition of fat is regulated by lncRNA. Fat has an important effect on meat quality. However, there are relatively few studies on lncRNAs in the subcutaneous adipose tissue of Duolang sheep and Small Tail Han sheep. In this study, RNA-Seq technology and bioinformatics methods were used to identify and analyze the lncRNA and mRNA in the subcutaneous adipose tissue of the two breeds of sheep. The results showed that 107 lnRNAs and 1329 mRNAs were differentially expressed. The differentially expressed genes and lncRNA target genes were significantly enriched in the biosynthesis of unsaturated fatty acids signaling pathway, fatty acid metabolism, adipocyte differentiation and other processes related to fat deposition. Among them, LOC105616076, LOC114118103, LOC105607837, LOC101116622, and LOC105603235 target FADS1, SCD, ELOVL6, HSD17B12 and HACD2, respectively. They play a key regulatory role in the biosynthesis of unsaturated fatty acids. This study lays a foundation for the study of the molecular mechanism of lncRNA on fat development, and has reference value for studying the differences in fat deposition between Duolang sheep and Small Tail Han sheep.

stages, and found the candidate lncRNAs regulate fat deposition by target genes. Cheng et al. 16 compared differences between preadipocytes and mature adipocytes by whole-transcriptome sequencing and constructed systematically regulatory networks according to the relationship predicted among the differentially expressed RNAs. They found that the lncRNAs play important roles in the regulatory networks and influence adipocyte differentiation. These studies indicate that lncRNAs might determine fat deposition and regulate adipogenic differentation. However, there are few studies on the expression profile and function of lncRNA in the subcutaneous adipose tissue of Duolang sheep and Small Tail Han sheep with differences in fat deposition.
Duolang Sheep 17 is a sheep breed with excellent meat-fat quality in Xinjiang, China. It has the characteristics of high meat yield and strong fecundity. Its meat is tender and juicy, with uniform fat deposition and no taint. Small Tail Han Sheep 18-20 is a dominant sheep breed in northern China, with high fecundity, early sexual maturity, stable genetic performance, and strong adaptability.Although the Small Tail Han sheep has many advantages, the meat body size is not obvious and the carcass meat production rate is low 21 . In addition, the tail of the Duolang sheep is large and the tail of Small Tail Han sheep is short and round. The two samples belonged to fat-tailed sheep. These two kinds of sheep provide good research objects for adipogenic differentiation and fat deposition. In this study, RNA-Seq technology and bioinformatics methods were used to sequence the transcriptome of lncRNA and mRNA in the subcutaneous adipose tissue of Duolang sheep and Small Tail Han sheep. In addition, we conducted a comprehensive analysis of them to have a deeper understanding of the role of lncRNA in fat deposition and lipid metabolism. These can lay the foundation for revealing the fat development of two breeds of sheep.

Materials and methods
Experimental animals and sample preparation. In this study, three healthy female Duolang sheep and three healthy female Small Tail Han sheep were fed a diet formulated to meet current nutritional requirements. All were 2 years old. All animals can drink and eat freely under natural light. The weights of the species were similar (50 kg), and all were healthy and in good physical condition. We collected subcutaneous adipose tissue samples located at the backfat. To ease the pain,we stunned them with electricity and then slaughtered. The subcutaneous adipose tissue was sampled into a 5 ml tube within 30 min after slaughter and immediately frozen in liquid nitrogen. Then transferred to refrigerator at -80° for long-term preservation and further total RNA extraction.
Isolation of total RNA and quality control. Total RNAs were extracted from the same amount of subcutaneous adipose tissue using TRIzol 22 (Invitrogen Life Technologies, Carlsbad, USA) according to the manufacturer's instructions. In addition, genomic DNA was removed using rDNAIRnase-free (TaKara). RNA quality was verified using 2100 Bioanalyzer (AgilentTechnologies, Santa Clara, CA, USA) and NanoDrop 2000 (Thermo Scientific, USA). Qualified RNA was used for sequencing library construction, and their OD260/280 ≥ 1.8, OD260/230 ≥ 1.0, RIN ≥ 8, the brightness of the 28S was significantly higher than the 18S, and the RNA concentrations of all samples were 200 ng/μL. cDNA library construction and RNA sequencing. The qualified total RNAs were used for library construction. RNA-seq transcriptome strand library was prepared following TruSeq stranded total RNA Kit from Illumina (San Diego, CA) using total RNA. The steps included removal of ribosomal (rRNA) and enrichment of mRNA. Then SuperScript double-stranded cDNA synthesis kit (Invitrogen, CA) was used and firststranded cDNA was synthesized with random hexamer primers. We removed the RNA template and synthesized a replacement strand, incorporating dUTP in place of dTTP to generate ds cDNA. This was followed by the addition of the end repair mix to blunt the sticky end, followed by the addition of an A base at the 3' end to form a Y-linker. The products after the adaptor were purified and sorted. The sorted products were used for PCR amplification and purification to obtain the final library. High throughput sequencing was conducted using the Illumina NovaSeq 6000 sequencing platform.
Reference genome mapping and transcriptome assembly. Clean data were obtained using fastp software (V0.19.5) filtering out the raw data containing joints, low-quality reads, N rate (N represents uncertain base information), higher sequences and too short sequences 23 , the remaining sequences were used for further analysis. Comparing clean reads with the reference genome GCF_002742125.1 using Hisat2 24 (V2.1.0). Analyzed valid data mapped reads. Then, we spliced mapped reads using StringTie 25 software (V1.3.3b). Comparing with the original genome annotation information, we found unannotated transcription regions and discovered new transcripts and new genes of the species.
Identification of potential lncRNA candidates. LncRNA is a non-coding RNA with a length of more than 200 nucleotides. Based on these features, we selected the transcripts with class symbols "x", "i", "j", "u", "o". On this basis, the lncRNAs with length ≥ 200 nt, number of exon ≥ 2 and ORF ≤ 300 bp were selected as candidates for preliminary screening. The absence of protein-coding ability was a key condition for judging lncRNAs. We used CPC 26 , CNCI 27 and Pfam 28 to eliminate transcripts with coding ability. Therefore, we obtained the newly predicted lncRNA in the fat sample, and identified known lncRNAs from NCBI and Ensembl databases.
Screening of differentially expressed genes and lncRNAs. The expression level of lncRNAs and mRNAs were corrected by the transcripts per million reads (TPM) method, and the gene length and sequencing depth were uniformized. Therefore, that the expression levels in different samples were consistent. The expres- www.nature.com/scientificreports/ sion levels between genes can be more intuitively. For experiments with biological replicates, use the DESeq2 based on negative binomial distribution to perform statistical analysis on raw counts, and obtain differential expression based on p value < 0.05, |log2FoldChange|> 1.
GO and KEGG enrichment analyse of differentially expressed genes. Gene Ontology Consortium (GO) can be used for functional enrichment analysis of differentially expressed genes, divided into cellular component (CC), molecular function (MF) and biological process (BP) 29,30 . Kyoto Encyclopedia of Genes and Genomes (KEGG) 31 can link genomic and functional information. Classification of functions was exercised 32 . We displayed the genes in the gene set on the KEGG pathway diagram, and display the KEGG annotation pathway diagram they participate in. The multiple test correction method Benjamini and Hochberg (BH) was used to correct the pvalue, resulting in padjust. When the padjust < 0.05, it was significantly enriched.
Protein-protein interaction (PPI) network analysis for differentially expressed genes. According to the STRING database (http:// strin gdb. org/), the PPI network analysis of differentially expressed genes was carried out. The relationship between the differentially expressed genes in the subcutaneous adipose tissue of Duolang sheep and Small Tail Han sheep was further studied. The visually analyze was conducted using Cytoscape software 33 with the obtained differential gene encoding protein interaction network data files, thereby obtaining key genes.
Prediction and functional analysis of the target genes of differentially expressed lncR-NAs. LncRNA is a non-coding RNA whose function is mainly to regulate target genes. According to the different mode of action of lncRNA, it is divided into homeopathic regulation (cis-regulate) and trans-regulation (trans-regulate). Homeopathic regulation and trans regulation of distant protein-coding genes were mentioned in Wang et al. 34 . The co-expression relationship between lncRNA and mRNA was realized by Pearson correlation coefficient (PCC) calculation. The setting standard was |PCC|> 0.8 and p value < 0.05 to screened the coexpressed lncRNA-mRNA. When screening target genes, genes within 100 kb of lncRNA were considered as target genes for cis 35 . Genes with |PCC|> 0.9 and a significant pvalue less than 0.05 in the co-expression analysis were used as target genes for the trans-regulate of differential lncRNAs. The potential role of lncRNA was studied through GO annotation and KEGG pathway enrichment analysis of target genes.

Real-time fluorescence-based quantitative PCR (qRT-PCR) verification. QRT-PCR was used
to verify the reliability of the sequencing results. Seven differentially expressed mRNAs and six differentially expressed lncRNAs were randomly selected. These genes come from PPI, lncRNA-mRNA network and fat metabolism related pathways. Three biological replicates were employed for each gene. 0.5 μg RNA was taken to synthesize cDNA template through GeneAmp PCR System 9700 (APPlied Biosystems, USA). The qRT-PCR analysis was conducted using LightCycler 480IIReal-time PCR Instrument (Roche, Swiss). The reaction system consisted of 1 μL cDNA, 5 μL of 2 × PerfectStart Green qPCR SuperMix, 0.2 μL of 10 μM forward primer, 0.2 μL of 10 μM reverse primer, and 3.6 μL of nuclease-free H2O. Reaction conditions were as follows: 94 °C for 30 s, 94 °C for 5 s, 60 °C for 30 s, 45 cycles. After the cycle, the melting curve was used to detect the specificity of the product: the temperature was slowly increased from 60 to 97 °C. The relative expression levels of genes between samples were calculated using 2−ΔΔCt method 36 . Data obtained were analyzed using GraphPad Prism (V8.0.1).
The student t-test (p < 0.05) was used for mean comparisons. All results were presented in bar charts with the means and their standard deviation (± SD). Primers used for the qRT-PCR are listed in supplementary table 6.

Statistical analysis.
All the data were presented as means ± SD. When comparisons were made, a Student's t-test was performed and p value < 0.05 was considered as statistically significant.
Ethics statement. All the procedures involving animals were approved by the animal care and use com-  (Fig. 1A). Most lncRNAs have 2 to 4 exons, while mRNAs have 6 or more exons than lncRNAs. Most lncRNAs have 19 exons (Fig. 1B). In Fig. 1C, the comparison shows that the lengths of lncRNA and mRNA have the same increasing and decreasing trend. However, the ratio of longer lncRNA is lower than that of mRNA. It is predicted that the open reading frame length of lncRNA is generally shorter than that of mRNA, as shown in Fig. 1D.
Expression level of genes and differential expression analysis. According to the transcript expression level measurement index TPM value, a violin plot was constructed, and the transcript expression levels in the two samples were analyzed separately (Fig. 1E) Table 2). In Fig. 2, the regulation of multicellular organismal process and the positive regulation of multicellular organismal process are the two most important enrichment processes in biological processes, followed by fatty acid derivative metabolic process, regulation of fat cell differentiation and significantly enrichment of fatty-acyl-CoA in metabolic process. In molecular functions, it is mainly enriched in items such as cargo receptor activity and carbohydrate binding. In cell components, it is mainly enriched in items such as membrane part and intrinsic component of membrane. Through GO enrichment analysis, it is found that the differentially expressed genes mainly regulate the signal transduction and lipid metabolism processes of the adipocytes of Duolang sheep and Small Tail Han sheep.
KEGG enrichment analysis of differentially expressed genes. The results of KEGG enrichment analysis showed that a total of 943 out of 1329 differentially expressed genes were annotated with KEGG. They were annotated into 306 signaling pathways, of which 20 pathways were significantly enriched (Supplementary Table 2). In Fig. 3, pathway enrichment results show that differentially expressed genes are enriched in a large number of pathways related to lipid metabolism, including biosynthesis of unsaturated fatty acids, arachidonic acid metabolism, glycerolipid metabolism, etc. In addition, differentially expressed genes are also significantly enriched in cell signal transduction pathways such as PI3K-Akt and TGF-beta signaling pathway to regulate proliferation, differentiation and apoptosis. In addition, the differentially expressed genes were also significantly enriched in insulin resistance related diseases caused by abnormal lipid metabolism. The above results indicate that the differentially expressed genes in the subcutaneous adipose tissue of Duolang sheep and Small Tail Han sheep are involved in multiple pathways related to lipid metabolism. Among them, genes enriched in biosynthesis of unsaturated fatty acids and arachidonic acid metabolism may be beneficial to Duolang sheep. It plays an important role in regulating the metabolism of subcutaneous fat in Small Tail Han sheep.

PPI network for differentially expressed genes. Fat deposition is the result of fat cell differentiation
and lipid metabolism. In this study, a PPI network of differentially expressed genes related to adipocyte differentiation and lipid metabolism was constructed through PPI network analysis and gene differential expression. In addition, genes that regulate Duolang sheep fat deposition and Small Tail Han sheep were determined (Supplementary Table 3). In the PPI network, combine-score ≥ 0.6 and degree ≥ 10 were used as thresholds to screen key nodes. Among them, SCD, DHCR24, PTGS2, TGFB1, FADS1 and other encoded proteins were at the key node positions (Fig. 4), which may be involved in subcutaneous fat metabolism. This has an important regulatory role in adipogenic differentiation.

Target genes of lncRNAs and functional analysis.
The differentially expressed lncRNAs in the subcutaneous adipose tissue of Duolang sheep and Small Tail Han sheep were used to predict the target genes of cis and trans. The results showed that 18 differentially expressed lncRNAs were predicted to target 28 genes in line with   Table 2). The prediction results of trans target genes showed that 56 of the 107 differentially expressed lncRNAs are co-expressed with differential mRNAs (Supplementary Table 4). Combining GO and KEGG analysis to screen out the genes related to adipocyte differentiation, lipid metabolism and related diseases, and screen with |PCC|≥ 0.9 as the threshold to construct a coexpression network for the trans-target genes of differential lncRNAs..We found that lncRNA may correspond to multiple mRNAs, and one mRNA may correspond to multiple lncRNAs, and the relationship between the two is not necessarily one-to-one. In Fig. 5, LOC101116622, LOC105603235, LOC114110986, LOC114114983, LOC114118103, LOC114108859, LOC114113946, LOC105614707, LOC105616344 are key lncRNAs. The function of the lncRNA can be inferred from the function of mRNA. The GO annotation and KEGG pathway enrichment of the target gene showed that it is significantly enriched in biological processes, including regulation of fat cell differentiation, lipid metabolism, fatty acid metabolism, and arachidonic acid metabolism, lipid metabolism pathway such as biosynthesis of unsaturated fatty acids. This indicates that the differentially expressed lncRNAs may regulate adipose tissue by participating in the above-mentioned signal pathways and biological processes. The target genes were selected to construct a lncRNA-mRNA network in biosynthesis of unsaturated fatty acids (Fig. 6) (Supplementary Table 5), and to explore the lncRNA-mRNA regulation mode.
qRT-PCR verification. QRT-PCR is considered as the golden standard for quantitative analysis of genes.
RNA-Seq data is significantly correlated with qRT-PCR results. Therefore, 13 differentially expressed genes were randomly selected to verify the RNA-Seq results. In Fig. 7, LOC101113583, LOC105614707, PPP2R5A, LOC114108859, PCK1, FADS1, LOC114117814 and PTGS2 were up-regulated in subcutaneous adipose tissue of Duolang sheep. COL1A1, LOC105605445, AKT3, LOC114116830 and LOC114112974 were up-regulated in subcutaneous adipose tissue of Duolang sheep. These results were consistent with the sequencing results, indicating the reliability of sequencing results.

Discussion
Fat deposition is regulated by many transcription factors, key genes and signaling pathways. This study applied RNA-Seq technology to analyze gene expression and reveal its biological characteristics. Duolang sheep and Small Tail Han sheep are high-quality local sheep breeds in China. In order to explore their differences in fat deposition, we use RNA-Seq technology to sequence the subcutaneous fat tissue. Then, we screened out the differentially expressed genes, and performed GO annotation and KEGG enrichment analysis on the 1329 differentially expressed genes. Finally, we screened out the target genes for fat deposit effects through analysis. In addition, the 107 differentially expressed lncRNAs and differentially expressed mRNAs were screened for co-expression analysis. Genes related to adipocyte differentiation, lipid metabolism and related diseases were screened out as the target genes of lncRNAs to reveal the function of lncRNAs. These genes will guide the breeding of high-quality livestock and poultry, improve meat quality, and provide potential targets for the treatment of fat metabolism-related diseases.   39,40 . In this study, unsaturated fatty acids biosynthesis was significantly enriched with 8 differentially expressed genes, including TECR, ACOT7, LOC114109406, FADS1, SCD, ELOVL6, HSD17B12, HACD2. These genes were regulated by 29 lncRNAs. TECR is a gene related to fatty acid synthesis, which would decrease the transcript level in NaF solution 41 . ACOT7 is an enzyme that catalyze the hydrolysis of acyl-CoAs to free fatty acids and CoA 42 . However, there were no report on TECR, ACOT7 and LOC114109406 gene regulating fat deposition in sheep. FADS1 and SCD are the key node position of PPI network. FADS1 is a member of the fatty acid dehydrogenase gene family and is related to all lipids 43 . It can encode Delta-5 (D5D) desaturase, which is the rate-limiting enzyme for the conversion of polyunsaturated fatty acids. This is considered as the main determinant of the level of polyunsaturated fatty acids 44 . That the expression of FADS1 in the lower adipose tissue of the cowhide is up-regulated under cold stimulation conditions. This promotes the production of polyunsaturated fatty acids in the lower adipose of the cowhide 45 . It was also confirmed in pigs and cattle 46,47 . In addition, FADS1 is also associated with metabolic diseases (obesity, metabolic syndrome) and cardiovascular diseases (arterial hypertension, coronary heart disease) 48,49 , In this study, FADS1 was up-regulated in Duolang sheep, which was consistent with previous study showing that FADS1 was highly expressed in subcutaneous adipose tissue of other species. Meanwhile, LOC105616076 might target FADS1 and involve in unsaturated fatty acids biosynthesis, thus regulating lipid metabolism.
Stearoyl-CoA desaturase (SCD) is an endoplasmic reticulum binding enzyme, belonging to the dehydrogenase family. It is the key enzyme that catalyzes the formation of unsaturated fatty acids from saturated fatty acids, especially palmitoleic acid and oleic acid. The high expression of SCD can promote the fat production, and the expression product can participate in the differentiation of pre-adipocytes and cell metabolism. There are many  Figure 5. Target relationship between lncRNA and mRNA. The circle in the node represents mRNA, and the triangle represents lncRNA, edge represents interaction between lncRNA and mRNA. In addition, red represents up-regulation of Duolang sheep, green represents down-regulation. www.nature.com/scientificreports/ subtypes of SCD, among which SCD1 is a downstream gene regulated by SREBP1. When the sterol regulatory element on the SCD1 promoter is combined with the SREBP1c transcription factor, transcription is upregulated and lipid synthesis increases 50 . The backfat thickness of Duolang was higher than that of Small Tail Han sheep. In this study, SCD is up-regulated in the subcutaneous fat tissue of Duolang sheep, which is consistent with previous studies in sheep 51 , cattle 52 and pig 53 . These results suggest that SCD promotes subcutaneous fat deposition in Duolang sheep and the mechanisms of SCD regulating fat deposition in cattle and pig might be same. In addition, LOC114118103 might participate in the biosynthesis of unsaturated fatty acids pathway by regulating SCD, and regulate the deposition of subcutaneous adipose tissue in Duolang sheep and Small Tail Han sheep. Long-chain fatty acid elongase 6 (ELOVL6) is one of the important regulatory factors that regulate fatty acid synthesis. It is the rate-limiting enzyme for the elongation of long-chain fatty acids. It mainly regulates cell metabolism, proliferation or apoptosis by affecting the composition of fatty acids and metabolites. The ELOVL6 polymorphism is significantly related to subcutaneous fat deposition in chickens 54 . The expression of the ELOVL6 requires activation of the retinoid X receptor RXR 55 . Its product stearic acid regulates mitochondrial function through transferrin receptor 1 (TFR1) 56 . Lin et al. 57 analyzed indicated that SCD and ELOVL6 expressions were positively correlated with the concentrations of polyunsaturated fatty acids in yak subcutaneous fat. However,  www.nature.com/scientificreports/ there is no report on ELOVL6 gene regulating fat deposition in sheep. In this study, ELOVL6 was used as the target gene of LOC105607837 to regulate fat deposition. HSD17B12 is a multifunctional isoenzyme, which plays an important role in the elongation of long-chain fatty acids. Increased gene transcription of HSD17B12 may lead to increased lipid synthesis. This can reduce 3-ketoacyl-CoA in the prolongation of VLCFA biosynthesis. It can be catalytically reduced to 3-hydroxyacyl-CoA 58,59 . In this paper, LOC101116622 regulated lipid metabolism by targeting HSD17B12. HACD2 is the main 3-hydroxyacyl-CoA dehydratase. It produces HACD2-deficient cells in the mammalian system, resulting in the reduction of ≥ C18 saturated and monounsaturated FAs 60 . This proves that HACD2 contributes to fat deposition. In this paper, LOC105603235 regulated the metabolism of fatty acids by targeting the HACD2, resulting in thicker backfat in Duolang sheep. The above studies, have shown that through the identification and analysis of genes in the subcutaneous adipose tissue of Duolang sheep and Small Tail Han sheep, a number of candidate genes involved in the biosynthesis of unsaturated fatty acids have been screened, including TECR, ACOT7, LOC114109406, FADS1, SCD, ELOVL6, HSD17B12 and HACD2. These genes may have reference value for studying the difference in fat deposition between Duolang sheep and Small Tail Han sheep.

Conclusions
In this study, RNA-Seq technology and bioinformatics methods were applied to explore the mechanism of fat deposition in sheep. We identified the differentially expressed lncRNAs and genes of subcutaneous adipose tissue between Duolang and Small Tail Han sheep. On this basis, 107 lncRNAs and 1329 mRNAs were differentially expressed. Combined with GO and KEGG analysis of differentially expressed genes, lncRNA participates in the biosynthesis of unsaturated fatty acids pathway through targeted mRNA. By this way, they regulated subcutaneous fat deposition in Duolang and Small Tail Han sheep. LOC105616076, LOC114118103, LOC105607837, LOC101116622 and LOC105603235, which might target FADS1, SCD, ELOVL6, HSD17B12 and HACD2, respectively, play key regulatory roles. Previous studies on lncRNAs have mostly appeared in model animals such as mice, or in pigs that are genetically similar to humans, and there is relatively few researches on lncRNAs in sheep. This study can provide useful information for understanding molecular mechanism of fat deposition in sheep. In addition, this study can help breed sheep with high meat quality as well as prevent and treat disease associated with fat metabolism.